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Abstract. These lecture notes discuss several topics in the physics of 
cosmic structure formation starting from the evolution of small-amplitude 
fluctuations in the radiation-dominated era. The topics include relativis- 
tic cosmological perturbation theory with the scalar-vector-tensor decom- 
position, the evolution of adiabatic and isocurvature initial fluctuations, 
microwave background anisotropy, spatial and angular power spectra, the 
cold dark matter linear transfer function, Press-Schechter theory, and a 
brief introduction to numerical simulation methods. 



1. Introduction 

The advent of precision cosmology with measurements of the Cosmic Microwave 
Background (CMB) radiation has brought a new focus on the physics of the z > 
1100 universe. The connection between two fossil relics of the early universe — 
CMB fluctuations and galaxies — offers the possibility for improved understand- 
ing of both as well as tighter constraints on cosmological parameters. These 
lectures present an overview of the physics of structure formation starting af- 
ter the inflationary epoch, going through recombination, and extending to the 
present day. 

Cosmological perturbation theory is often regarded as a highly technical 
subject beyond the scope of graduate courses in cosmology. However, these lec- 
tures may suggest how, with minor simplification, it can be integrated into a 
treatment of large scale structure. The potential benefits are many — a detailed 
understanding of CMB anisotropy, a quantitative explanation for the key fea- 
tures of the cold dark matter transfer function, and an extension of the intuition 
developed from the Newtonian theory of large scale structure to the physics of 
the post-inflationary universe. The aim of the current notes is more modest, 
but it is hoped that they provide some direction for those wishing to delve more 
deeply into cosmological perturbation theory and structure formation. 



2. Overview of Relativistic Cosmological Perturbation Theory 

The goal of cosmological perturbation theory is to relate the physics of the early 
universe (e.g. inflation) to CMB anisotropy and large-scale structure and to 
provide the initial conditions for numerical simulations of structure formation. 
The physics during the period from the end of inflation to the beginning of 
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nonlinear gravitational collapse is complicated by relativistic effects but greatly 
simplified by the small amplitude of perturbations. Thus, an essentially complete 
and accurate treatment of relativistic perturbation evolution is possible, at least 
in the context of simple fluctuation models like inflation. Numerous articles 
have been written on this topic. See [l], for references and for more detailed 
treatments similar in style to the present lecture notes. 

The starting point for cosmological perturbation theory is the metric of a 
perturbed Robertson- Walker spacetime, 



ds 2 



o 2 (t) -dr 2 + r y ij (x)dx l dx j + h^(x,T)dx tl dx u . (1) 



Spatial coordinates take the range 1 < i,j < 3; x % (or x for all three) is a 
comoving spatial coordinate; r is conformal time; a(-r) is the cosmic expansion 
scale factor; units are chosen so that the speed of light is unity; and "fij{x ) is the 
3-metric of a maximally symmetric constant curvature space. Conformal time 
is related to the proper time measured by a comoving observer (i.e. one at fixed 
x) by dt = a{r)dr. For cosmological perturbation theory it is more convenient 

than proper time. The metric perturbations are given by h^ u = g$ /a 2 . 

The Robertson- Walker model is characterized by one function of time, a(r), 
and one constant, the curvature constant K. The expansion factor obeys the 
Friedmann equation, 

h -) 2 = ^Ga 2 p{a)-K (2) 
a J 3 

where a dot denotes d/dr and p(a) is the total mean energy density. Note that 
the usual Hubble parameter is H = a /a 2 . 

The spatial part of the Robertson- Walker metric takes a simple form in 
spherical coordinates (x,0,(p): 

j ij dx i dx j = d X 2 + r 2 { X ) (d6 2 + sin 2 ## 2 ) . (3) 

The angular radius r(x) takes one of three forms depending on the sign of K: 

(-t=sm( X VK), K>0, 
r(x) = {x, K = 0, (4) 

[-^sinhix^K), K<0. 

These cases are referred to as closed, flat, and open, respectively. Many authors 
choose units of length such that x = il or 0- However, it is convenient to let K 
have units (inverse length squared). From the Friedmann equation, K = (Qq — 
1)Hq where the subscript denotes the present value and is the total density 
parameter (including dark energy). Note that the range of x is [0>°°) m t ne 
flat and open models but is [0,ir/y/~K\ in the closed model. The flat Robertson- 
Walker model favored by inflation is conformal to Minkowski spacetime, i.e. the 
metrics are identical up to an overall conformal factor a 2 (7~). 
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2.1. Scalar- Vector- Tensor Decomposition 

In linear perturbation theory, the metric perturbations L v are regarded as a 
tensor field residing on the background Robertson- Walker spacetime. As a sym- 
metric 4x4 matrix, h^ u has 10 degrees of freedom. Because of the ability to 
make continuous deformations of the coordinates, 4 of these degrees of freedom 
are gauge- (coordinate-) dependent, leaving 6 physical degrees of freedom. A 
proper treatment of cosmological perturbation theory requires clear separation 
between physical and gauge degrees of freedom. 

The metric degrees of freedom in linear pertubation theory were classified 
by Lifshitz in 1946 Q. Lifshitz presented the scalar- vector-tensor decomposition 
of the metric. It is based on a 3 + 1 split of the components, which we rewrite 
as follows: 

/i o = -2V> , h 0i = Wi , hij = 2 {(jrfij + Sij) with 7^ £„■ = 0. (5) 

Here, 7 y is the matrix inverse of 7^. The trace part of hij has been absorbed 
into 4> so that Sij has only 5 independent components. 

With the space-time (3 + 1) split, we use 7^' (or 7^) to lower (or raise) 
indices of spatial 3-vectors and tensors. For convenience, spatial derivatives will 
be written using the 3-dimensional covariant derivative Vi defined with respect 
to the 3-metric 7™; it is a three-dimensional version of the full 4-dimensional 
covariant derivative presented in general relativity textbooks. For example, if 
K = 0, we can choose Cartesian coordinates such that 7^' = 5ij and Vj = d/dx l . 

The scalar-tensor-vector split is based on the decomposition of a vector 
into longitudinal and transverse parts. For any three- vector field Wi(x), we may 
write 

Wi = wj + wf~ where V x w " = V • w ± = . (6) 

The curl and divergence are defined using the spatial covariant derivative, e.g. 
V • w = j^ViWj. 

The longitudinal/transverse decomposition is not unique (e.g. one may 

always add a constant to w\ ) but it always exists. The terminology arises because 

in the Fourier domain w\ is parallel to the wavevector while wf- is transverse 

(perpendicular to the wavevector). Note that t«J = Vi<p w for some scalar field 
(ft w . Thus, the longitudinal/transverse decomposition allows us to write a vector 
field in terms of a scalar (the longitudinal or irrotational part) and a part that 
cannot be obtained from a scalar (the transverse or rotational part). 

A similar decomposition holds for a two-index tensor, but now each index 
can be either longitudinal or transverse. For a symmetric tensor, there are 
three possibilities: both indices are longitudinal, one is transverse, or two are 
transverse. These are written as follows: 

= &ij + $>ij + Sjj i (7) 

where 

l 3k VkS lj = 7 jk V k sl + ^ k V k S^j . (8) 

The first term in equation (||) is a longitudinal vector while the second term is 
a transverse vector. The divergence of the doubly-transverse part, S^, is zero. 
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For a traceless symmetric tensor, the doubly and singly longitudinal parts can 
be obtained from the gradients of a scalar and a transverse vector, respectively: 

4 = - \^ 2 ) <t>s, S± = V t Sf + VjSi- . (9) 

Although we will not prove it here (see [IJ for the details), under an in- 
finitesimal coordinate transformation sh and Sh can change while S?j is in- 
variant. Similarly, the longitudinal and transverse parts of W{ = hoi are both 
gauge-dependent. 

Now we have the mathematical background needed to perform the scalar- 
tensor- vector split of the physical degrees of freedom of the metric. The "tensor 
mode" represents the part of hij that cannot be obtained from the gradients of 
a scalar or vector, namely Sjj. The tensor mode is gauge-invariant and has two 
degrees of freedom (five for a symmetric traceless 3x3 matrix, less three from 
the condition 7 jfc VfcS^ = 0). Physically it represents gravitational radiation; 
the two degrees of freedom correspond to the two polarizations of gravitational 
radiation. Gravitational radiation is transverse: a wave propagating in the z- 
direction can have nonzero components h xx — h yy and h xy = h yx but no others. 
The tensor mode SJa behaves like a spin-2 field under spatial rotations. Note 
that my definition of the tensor mode strain Sij here differs by a factor of 2 from 
the conventional hj^ in Q. 

The "vector mode" behaves like a spin-1 field under spatial rotations. It 
corresponds to the transverse vector parts of the metric, which are found in wf~ 
and S^j. Each part has two degrees of freedom. Although we will not prove it 
here (see [Q] for the details), it is possible to eliminate two of these degrees of 
freedom by imposing gauge conditions (coordinate conditions) . The synchronous 

gauge of Lifshitz |J is one popular choice, with wf- = wj = 0. Another choice 
is the "Poisson" or transverse gauge ]|] which sets S^- = 0. In either case, 
there are two physical degrees of freedom and they correspond physically to 
gravitomagnetism. Although this effect is less well known than gravitational 
radiation or Newtonian gravity, it produces magnetic-like effects on moving and 
spinning masses, such as the precession of a gyroscope in the gravitational field 
of a spinning mass (Lense-Thirring precession). This phenomenon has not yet 
been discovered experimentally but should be measured by the Gravity Probe 
B satellite to be launched in 2002 g. 

The "scalar mode" is spin-0 under spatial rotations and corresponds physi- 
cally to Newtonian gravitation with relativistic modifications. The scalar parts 

of the metric are given by (/), tp, w\ , and S-j. Any two of these may be set to zero 

by means of a gauge transformation. One popular choice is w\ = sfj = 0, also 
known as the conformal Newtonian gauge. It corresponds to the scalar mode in 
the transverse gauge, defined by the gauge conditions 

^V tWj = , ^ k V k S tj = . (10) 

This gauge in linearized general relativity is the gravitational analogue of Coulomb 
gauge in electromagnetism, where the magnetic vector potential is transverse, 
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ViA 1 = 0. In the gravitational transverse gauge, Wi is transverse and Sij is dou- 
bly transverse. This gauge is convenient for developing intuition although not 
necessarily the best for computation. The variables ((/), tp, w^, S]j) correspond 
to gauge-invariant variables introduced by Bardeen || as linear combinations of 
metric variables in other gauges. 

In linear perturbation theory, the scalar, vector, and tensor modes evolve 
independently. The vector and tensor modes produce no density perturbations 
and therefore are unimportant for structure formation, although they do perturb 
the microwave background. In the remainder of these lectures, only the scalar 
mode will be considered. 

The Einstein equations give the equations of motion for the metric per- 
turbations in terms of the energy-momentum tensor, the source of relativistic 
gravity. Here we consider the case of a perfect fluid (or several perfect fluid 
components combined), for which 

T*" = (p + p)7"y + pg» v (11) 

where p and p are the proper energy density and pressure in the fluid rest frame 
and is the fluid 4-velocity. 

We split T^ u into time and space components as we did for the metric. It 
is convenient to use mixed components because the conformal factor a 2 then 
cancels, yielding 

T° = -p(x,r) = -[p(r) + 5p(x,T)] , 

A = [p(r)+p(r)}v l (x,T) = -[p(T)+p(T)]V t W , (12) 
T 3 = [p(t) + ^(£,t)]<5%. , 

where p(r) and p{r) are respectively the energy (or mass) density and pres- 
sure of the Robertson- Walker background spacetime, Vi = ^ijdx 3 /dr is the fluid 
3- velocity (assumed nonrelativistic), and W is a velocity potential. We are as- 
suming that the velocity field is irrotational or, if it is not, that the matter 
fluctuations are linear so that the longitudinal and transverse velocity compo- 
nents evolve independently. This is a good approximation prior to the epoch of 
galaxy formation. 

We are also making another approximation, namely that the shear stress 
(the non-diagonal part of T l •) is negligible compared with the pressure. This 
is a bad approximation for relativistic neutrinos, whose free-streaming leads to 
a non-isotropic momentum distribution after neutrino decoupling at a tempera- 
ture of about 1 MeV. (A similar but smaller effect occurs for photons after they 
decouple at a temperature of about 0.3 eV.) However, the effect of this error 
is modest during the radiation-dominated era (photons contribute more pres- 
sure than neutrinos, and are isotropic) and is unimportant during the matter- 
dominated era, when the nonrelativistic matter density is the major source of 
gravity. We introduce this approximation here in order to simplify the treat- 
ment; full calculations include the effects of neutrino shear stress Q. 

When the stress tensor is isotropic, the Einstein equations in transverse 
gauge imply that the two scalar potentials <\> and ip are equal Q and the metric 
takes the simple form 

ds 2 = a 2 {T) + 2(f))dT 2 + (1 - 2(f))-ii j dx l dx j \ . (13) 
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This is a simple cosmological generalization of the standard weak-field limit of 
general relativity, which is recovered by setting a = 1 and 7« = Sy. 

Given the stress-energy tensor, we can now obtain field equations for 4>(x, r) 
using the Einstein equations. The only independent equation for the Robertson- 
Walker background is the Friedmann equation (0). For the perturbations, we 
obtain 



(V 2 + 3K)(f> = AirGa 2 



a 



5p + 3-(p + p)W 



a 



(14) 



d T (/) + -(/) = AnGa 2 (p + p)W , (15) 



a 



and 



dU + 3-d T (j) - (%nGa 2 p + 2K](b = 4TiGa 2 5p . (16) 
a v / 

There are more equations than in Newtonian gravity because the Einstein equa- 
tions have local energy-momentum conservation built into them. 

Equation ( |i~4|) is remarkably like the Poisson equation of Newtonian grav- 
ity, especially when one realizes that the factor of a 2 is needed to convert the 
Laplacian from comoving to proper spatial coordinates. Otherwise there are 
three differences. First, the spatial curvature adds a scale to the Laplacian in a 
curved space. Second, the density source is 5p = p — p rather than p because 
(ft represents a metric perturbation on top of the background Robertson- Walker 
spacetime. Third, momentum density (through W) is a source for gravity in 
general relativity. However, this term is comparable with the density term only 
on length scales comparable to or larger than the Hubble length H . 

Two of our field equations are unfamiliar from Newtonian dynamics. Equa- 



tion (15) is reminiscent of the Zel'dovich approximation in that the velocity 
potential is proportional to the gravitational potential in the linear regime when 
(f>(x,T) = <f>o(x )D ( / ) (t). However, it is a more general result whose Newtonian 
limit may be derived by combining the continuity equation with the time deriva- 
tive of the Poisson equation. 

Equation ( |i~6| ) is more surprising, as it gives the gravitational potential 
through a second-order time evolution equation as opposed to the action-at-a- 
distance Poisson equation (|14|). The reader may well wonder about the com- 
patibility of the two equations as well as the physical validity of the Poisson 



equation in general relativity, but equations (14)-(16) are all correct. See ]lj] 
for a full discussion. In brief, causality is restored by the tensor mode. The 
time evolution equation for the potential can be regarded as arising from local 
energy-momentum conservation (which is built into the Einstein equations) com- 
bined with the Poisson equation. Indeed, an equation similar to equation (16) 



can be derived in Newtonian gravity from the fluid and Poisson equations. Its 
utility arises from the fact that when pressure perturbations are gravitationally 
unimportant (i.e. when matter and/or dark energy dominate over radiation), 
the time evolution of perturbations can be obtained from a single differential 
equation depending only on the background curvature, pressure, and density. 

2.2. Perfect Fluid Model 

To conclude this section, we present a simple model for understanding CMB 
anisotropy and large scale structure. We reduce the matter and energy content 
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of the universe to a vacuum energy (with p v = A/8ttG contributing to the mean 
expansion rate but not to the fluctuations) and two fluctuating components, 
CDM and radiation (photons and neutrinos, the former coupled to baryons by 
electron scattering until recombination). The CDM and radiation are treated 
as perfect fluids, ignoring the free-streaming of neutrinos and the diffusion and 
streaming of photons during and after recombination. 

The fluid equations follow from local energy-momentum conservation, V^T^ 
0. For a perfect fluid, with the scalar mode only (no vorticity and no gravita- 
tional radiation), the fluid equations in the metric of equation (13) are 



d T p + 3(--d T 0){ P + p) + v i (p + pW 



o 



117) 



and 



9 T (P + P)V* 



+ 4-(p + p)v l + V iP +(p + p)V< 
a 



(18) 



These are familiar from Newtonian fluid dynamics, with some extra terms. The 
damping (a/a) terms arise from Hubble expansion and the use of comoving 
coordinates. The d T (f) term in the continuity equation ([17]) arises from the de- 
formation of the spatial coordinates, i.e. the <f>8ij contribution to the metric. 
(The flux term in eq. 17 describes transport relative to the comoving coordi- 
nate grid, which is deforming when d T <f) ^ 0.) The energy flux (or momentum 
density, they are equal) includes pressure in relativity because of the pdV work 
done by compression. Equations ([[7]) and ( |i~8| ) apply separately to perfect fluid 
component. 



3. Adiabatic and Isocurvature Modes 

This section discusses the evolution of perturbations from the end of inflation 
(taken to be at r = for all practical purposes) through recombination using a 
two-fluid approximation. We approximate the matter and energy content in the 
universe as being two perfect fluids: Cold Dark Matter (a zero-temperature gas) 
and radiation (photons coupled to electrons and baryons before recombination, 
plus neutrinos which are approximated as behaving like photons). The net 
energy density perturbation is thus 

5p(x, t) = p c (a)S c (x, t) + p r (a)5 r (x, r) , (19) 

where p c oc a -3 and, if we neglect the contribution of baryons, p r oc a -4 . We 
could add a third component for the baryons, but this complicates the presen- 
tation without adding essential new behavior. In the interest of economy, we 
leave baryons out. We define the relative density contrasts 5 C = 5p c /p c and 
5 r = Spr/ p r . 

The total pressure is simply p(x,r) = |p r (a)(l + 5 r ) + p v where p v = 
—p v is the spatially constant negative pressure of vacuum energy (cosmological 
constant), should any be present. 

As we noted in the previous section, density perturbations couple gravita- 
tionally only to the scalar mode of metric perturbations, and the scalar mode 
cannot generate transverse vector fields. As a result, the peculiar velocity fields 
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of CDM and radiation are longitudinal and are fully characterized by their po- 
tentials W c and W r for CDM and radiation, respectively. 

As a final simplifying assumption, we will assume that we are studying ef- 
fects on distance scales less than the curvature length |i^| -1 / 2 = | Oo — 1 
which is at least 5 Gpc and possibly infinitely large. 

With these assumptions, we can linearize the fluid equations (|l 
for 5 2 . <C 1 and 5 2 . <C 1 to obtain 



and ph 



d T 5 c = V 2 W C + 3d T c 



d T W c + -W c 
a 



d T s r = -V 2 W r + Ad T <p , 



d T W r 



:6r + 



(20) 



These equations are easy to understand. The d T (p terms have the same ori- 
gin as in equation (17), namely the deformation of the comoving coordinate 
system. The velocity potential terms are related to the velocity divergence, 
V ■ v = — V 2 W. The radiation density fluctuations grow at 4/3 the rate of 
matter fluctuations to compensate for the fact that the photon number density 

is proportional to pV^ ■ As expected, the velocity equations show that matter 
velocities redshift (for <j> = 0) with v c oc a -1 while radiation velocities (e.g. the 
speed of light) do not redshift. (Momenta redshift oc a -1 but the speed of light 
is unaffected by expansion.) Finally, radiation pressure gives rise to an acoustic 
restoring force — Vp r = — (l/3)/? r V5 r which shows up as the (l/4)5 r term in the 
radiation longitudinal velocity equation. 

Equations (p0|) are coupled by gravity, which can be determined from the 



Poisson equation (|14j), which becomes 



V 2 d> = 4vrGa 2 



(p r 5 r + PcS c ) + - (4:p r W r + 3p c W c 
a 



(21) 



We can solve equations (^) and (^T|) most easily by expanding the spatial 
dependence in eigenfunctions of the Laplacian V 2 . In flat space these are plane 
waves exp(ik ■ x), and this is a good approximation even if the background is 
curved, provided k 3> |i^| 1//2 . This approximation is valid for almost all appli- 
cations except large angular scale CMB anisotropy. Note that k is a comoving 
wavevector; the physical wavelength is 2ira/k. 

By expanding the spatial dependence in plane waves, we now have a fourth- 
order linear system of ordinary differential equations in time. We are free to 
define linear combinations of the variables (5 c ,5 r ,W c ,W r ) as our fundamental 
variables. In terms of the mechanisms for generating primeval perturbations, the 
most natural variables are the metric perturbation and the specific entropy 



V = 



5p — c 2 5p 3 
Pec 2 4 



(22) 



where 



c 2 = 



dp/dr 
dp/dr 



3 1 + 



3pc 

4p T 



(23) 
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is the effective one-fluid sound speed of the matter and radiation. (Acoustic 
signals actually do not propagate at this speed; they propagate with speed 3 -1 / 2 
through the radiation and they do not propagate at all through CDM.) 

For a two-component radiation plus CDM universe, the solution of the 
Friedmann equation (^) gives 

p c _ o(r) _ r / r \ 2 _ _ (a COL \ 1 / 2 x _ 19Mpc 



Pr ^eq T t 



The radiation-dominated era ends and matter-dominated era begins at redshift 
l + z eq = a~q = 2.5 x 10 4 Sl m /i 2 . A cosmological constant has no significant effect 
on y(r) provided that, during the times of interest, the vacuum energy density 
is much less than the radiation or cold dark matter density. 

The fluid and Poisson equations can now be combined to give a pair of 
second-order ordinary differential equations, 

V 




2yr| ' 

^ T n+-drV+^V = \y 2 k^. (25) 
ocj a 4 b 

It is interesting to note that (j> and r\ evolve independently aside from the source 
term each provides to the other. The coupling implies that entropy perturbations 
are a source for the growth of gravitational perturbations and vice versa. 

For a given wavenumber k, there are two key times in the evolution of cf> 
and rj: the sound-crossing time r ~ 7r/(£;c s ) and the time of matter-radiation 
equality at r w r e . [To be precise, a = a eq at r/r e = 2(V2 - 1).] For r < r e , <f> 
and rj decouple and they both have solutions that decay rapidly with r as well as 
solutions that are finite as r — ► 0. The latter are conventionally called "growing 
modes" while the former are called decaying. The growing mode solution in the 
radiation era r <C r e is 

^3 
4>(k,r) = - — ( sin ujt — lot cos lot) A(k) 
[lot) 6 



T 

+ — 



1 , ^ 
+ —, 777 - 7 7T (COS LOT + LOT Sin LOT) 



(cjt) 4 2(wt 



/(fc) ,(26) 



ri(k,T) = I(k) + 9 



In lot + C — Ci(u>r) + - (cos lot — 1) 



A(k) , 



where lo = k/\^3 is the phase of acoustic waves in the radiation fluid, C = 
0.5772 ... is the Euler constant, and Ci(x) is the cosine integral defined by 

Ci(a) = C + Inx + f X C ° St ~ 1 dt . (27) 
Jo t 

In equation (]26| ) we have neglected all terms 0(r/r e ) aside from the lowest-order 
contribution of I(k ) to 4>. 
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The general growing-mode solution contains two /c-dependent integration 
constants, A(k) and I(k), which represent the initial (r — > 0) amplitude of 
4> and 7], respectively. Conventionally they are called the "adiabatic" (so-called 
because the entropy perturbation 77 = at r = 0, although "isentropic" would be 
more appropriate) and "iso curvature" (so-called because <f> = at r = 0) mode 
amplitudes. In addition, there are decaying modes whose amplitude at r = r e 
is negligible unless the initial conditions were fine-tuned or there are external 
sources such as topological defects to continually regenerate the decaying mode. 
The solution of equations ( ]26| ) is valid only for inflation or other mechanisms that 
generate fluctuations exclusively in the early universe, but not for topological 
defects. 

As a result of reheating, inflation produces fluctuations with spatially con- 
stant ratios of the number densities of all particle species. For CDM and ra- 
diation, for example, p T ' / Pc = constant hence 1 = and thus 77 = until 
radiation pressure forces begin to separate the radiation from the CDM, which 
takes a sound-crossing time or r ~ u^ 1 . Thus, although inflation produces the 
adiabatic mode of initial fluctuations, the differing equations of state p(p) for 
the two-fluid system lead to entropy perturbations proportional to the initial 
adiabatic mode amplitude A(k ) (see the second of eqs. p6|). 

Physically, the number density of particles of all species are proportional to 
each other in the adiabatic mode as long as there has been too little time for 
particles to travel a significant difference via their thermal speeds. However, the 
differing thermal speeds of different particle species (e.g. photons and cold dark 
matter) eventually cause the particle densities to evolve differently, so that the 
number density ratios are no longer constant and the entropy perturbation is no 
longer zero. 

First-order phase transitions, on the other hand, change the equation of 
state without moving matter or energy on macroscopic scales, hence 5p r + 5p c = 
C = r = initially implying = initially but 77 7^ 0. This isocurvature mode 
is characterized by the initial entropy perturbation amplitude I(k ). As the first 
of equations (f26|) shows, the isocurvature mode generates curvature fluctuations 
when the universe becomes matter-dominated. 

The isocurvature CDM inflationary model has a fine balance such that 
o~r + y$c = initially. This may involve nonlinear fluctuations of 5 C at early 
times, but 6 r is negligible early in the radiation-dominated era y <C 1. For this 
reason, the isocurvature mode was originally called isothermal []|]. 

On small scales, kr S> 1, equations fl25| ) may be solved using a WKB 
approximation. For our purposes it suffices to take the ujt 3> 1 limit of equations 
©: 

COSWT . 1 f UJT ""\ , ~r* k , . 

0~-3^^2Mk) + — 2 --sinwr /Jj . 28 

The radiation density fluctuation scales as 

fc 2 2 2 
Or ~ ~ 2 - ~ k T <t> > 



which oscillates in proportion to cosujt for the adiabatic mode and grows in 
proportion to y ~ r/r e for the isocurvature mode. 
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Note from equation (^) that the isocurvature mode induces oscillations in 
the potential with a tt/2 phase shift relative to the adiabatic mode. This shift 
translates into a shift in the positions of the acoustic peaks in the CMB angular 
power spectrum. 

3.1. The Adbiabatic Growing Mode as a Standing Wave 

It is instructive to consider the general adiabatic perturbation in the radiation- 
dominated era including the decaying mode. This problem is easy to analyze 
and shows the role of acoustic waves in the radiation-dominated era. 

The two linearly independent solutions of the first of equations ( |25|) for 
y 1 and adiabatic initial conditions (rj = 0) are 

3 3 

7 (sin lot — lot cos ujt) , 4>d = ~, — ^(coswr + wT sinur) (29) 



where the subscripts label the growing and decaying modes, respectively. They 
are the real and imaginary parts of 

4> u = <f>a + %<f>d = T^yj(i ~ 0JT)e- iwT . (30) 

For a given spatial harmonic, exp(ik-x )<f> w represents a wave traveling in the 
fe-direction. Thus, the growing mode 4>g = 5 represents a superposition 

of waves traveling along k and — k, i.e. a standing wave 



b g cos kx = TZe 



2{0JTf 



(jT)e i{k -^ T) + H " iOT)e i{ - k - i+ulT) . (31) 



The waves travel at the sound speed u/k = l/\/3 in the radiation fluid. 

The entropy perturbation induced by a traveling wave potential perturba- 
tion is easily calculated in the radiation-dominated era using the Green's function 
method applied to the second of equations (^) with = <j> w for the source term. 
The result is 



9fl(l-e-)dx- 9 -(l-e—) . (32) 



The real part of this equation gives the A{k ) contribution to 77 in equations (|26|). 



4. CMB Anisotropy 

This section presents a simplified treatment of CMB anisotropy, with the aim of 
highlighting the essential physics without getting lost in all of the details. More 
complete treatments are found in 10]. 



The microwave background radiation is fully described by a set of photon 
phase space distribution functions. Ignoring polarization (a few percent effect), 
all the information is included in the intensity or in f(x,p,r) where fd 3 p is 
the number density of photons of conjugate momentum p at position and time 
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(x,t). The conjugate momentum is related to the proper momentum measured 
by a comoving observer, p/a, so that p = constant along a photon trajectory in 
the absence of metric perturbations. 

Remarkably, despite metric perturbations and scattering with free electrons, 
the CMB photon phase space distribution remains blackbody (Planckian) to 
exquisite precision: 



f(x,p,T) = /pianck (lJf< ) -A MancU 



P 



kT (l + A) 



(33) 



where To = 2.728 K is the current CMB temperature and A(x, n, r) is the 
temperature variation at position (x, r) for photons traveling in direction n. 
The phase space density is blackbody but the temperature depends on photon 
direction as a result of scattering and gravitational processes occurring along 
the line of sight. 

The phase space density may be calculated from initial conditions in the 
early universe through the Boltzmann equation 



Df _ df df dx df dp 
Dt dr dx dr dp dr 



df_ 

d,T 



(34) 



where the right-hand side is a collision integral coming from nonrelativistic 
electron-photon elastic scattering. In the absence of scattering, the phase space 
density is conserved along the trajectories 



dx _ p 
— = n=- 
dr p 



1 dp 
p dr 



-n 



V0 + d T (f> , — 



dn 
ah 



V - n(n- V) <f> . (35) 



Note that the photon energy changes via gravitational redshift (the first term in 
dp/dr) or a time-changing potential, but in either case the energy of all photons 
is changed by the same factor. The photon direction changes by transverse 
deflections due to gravity; dn/dr is twice the Newtonian value, a result well- 
known in gravitational lens theory [pL X[J - Because df /dn = for the unperturbed 
background (the anisotropy arises due to the metric perturbations), gravitational 
lensing affects the CMB anisotropy only through nonlinear effects on small scales 

111- 

The procedure for computing CMB anisotropy is to linearize the Boltzmann 
equation assuming A 2 <C 1 [13|. Until recently, the temperature anisotropy A 
was expanded in spherical harmonics and the Boltzmann equation was solved 
as a hierarchy of coupled equations for the various angular moments In 
1996, Seljak and Zaldarriaga Q developed a much faster integration method 
call CMBFAST based on integrating the linearized Boltzmann equation along 
the observer's line of sight before the angular expansion is made: 



A(n, r ) 



dxe 



-t(x) 



1 



n ■ V<j> + d T 4> + an e <JT ( -<5 7 + v e • n + pol. term I 



(36) 

where tq is the present conformal time, x is the radial comoving coordinate of 
equation (||), subscript "ret" denotes evaluation using retarded time r = To — x> 
(L is the relative density fluctuation in the photon gas (<L = 6 r in our two-fluid 
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approximation), v e is the mean electron velocity (i.e. the baryon velocity, which 
equals v r in our two-fluid approximation), <tt is the Thomson cross section, and 
the Thomson optical depth is 

rx 

tt(x) = / dx(an e a T ) ret . (37) 
J o 

We have left out small terms ( "pol. term" ) due to polarization and the anisotropy 
of Thomson scattering in equation j36l). 

The terms in equation ( |36| ) are easy to understand. The exp(— tt) factor 
accounts for the opacity of electron scattering, which prevents us from seeing 
much beyond a redshift of 1100, where tt ~ 1. The CMB anisotropy appears to 
come from a thin layer called the photosphere, just like the radiation from the 
surface of a star. 

The two gravity terms give the effective emissivity due to the photon energy 
change caused by a varying gravitational potential. For a blackbody distribution, 
if all photon energies are increased by a factor 1 + e, the distribution remains 
blackbody with a temperature increased by the same factor. Thus, the line- 
of-sight integration of the fractional energy change of equation (^) translates 
directly into a temperature variation. 

The terms proportional to the Thomson opacity an e o"T are the effective 
emissivity due to Thomson scattering. Bearing in mind that the photon-baryon 
plasma is in nearly perfect thermal equilibrium (the temperature variations are 
only about 1 part in 10 5 ), the photons we see scattered into the line of sight 
from a given fluid element have the blackbody distribution corresponding to the 
temperature of that element. Recalling that the energy density of blackbody 
radiation is p 7 oc T 4 , we see that is simply the fractional temperature 
variation of the fluid element. In other words, if we carried a thermometer 
to the photosphere and measured its reading in direction — n, it would read a 
fraction \8-y different from the mean. Now the energy density is defined in the 
fluid rest frame, while the fluid is moving with 3- velocity v e , so the temperature 
measured by a stationary thermometer (with fixed x l ) is changed by a Doppler 
correction v e ■ n. (If the photosphere is approaching us, v e \\ n, the radiation is 
blue-shifted.) 

Comparing the CMB with a star, when we look at the surface of a star we see 
the temperature of the emitting gas (to the extent that local thermodynamic 
equilibrium applies), corresponding to the term. For ordinary stars the 
Doppler boosting and gravitational redshift effects are negligible, although they 
are appreciable for supernovae and white dwarfs, respectively. For the CMB, 
on the other hand, all four emission terms in equation ( |36| ) are comparable in 
importance. 

It is instructive to approximate the photosphere as an infinitely thin layer 
by adopting the instantaneous recombination approximation, according to which 
the free electron fraction and hence opacity drops suddenly at r rec , the conformal 
time of recombination at z ~ 1100: 

r T ( X ) = fj°' X>X r ec = r -r rec , 

Substituting into equation ( |36|) and integrating by parts the gravitational red- 
shift term gives the result first obtained in another form by Sachs and Wolfe in 
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A(n) « Qj 7 + + 3 X W 7 ) + jf*"" d X (25 T 0) ret (39) 
where we hcLVG dropped, the argurnent To from A ; Xrcc — — ^rec? 

and W 7 is 

the velocity potential for the photons (W 7 = W r in the two-fluid approximation 
discussed above). We could have used the last of equations (|20| ) to write the 
sum of the intrinsic (j5<y) and gravitational ((f)) contributions as d T W r . The 
combination of derivatives (d T + d x )W r is the time derivative along the past 
light cone; the results are evaluated at recombination. 

While it is interesting that the primary anisotropy depends on the velocity 
potential derivative along the line of sight, there appears to be no fundamen- 
tal significance to the simple dependence in equation (|39|), aside from the fact 
that CMB anisotropy is produced by departures from hydrostatic equilibrium. 
(In hydrostatic equilibrium, v r = 0.) If the radiation gas were in hydrostatic 
equilibrium in a static gravitational potential, there would be no primary CMB 
anisotropy. Indeed, it can be shown from purely thermodynamic arguments that 
T(x ) oc exp[— cj)(x )] (hence |<5 7 + <fi = constant) for a photon gas in equilibrium. 

Sachs and Wolfe |14| showed that for adiabatic perturbations in the matter- 
dominated era, on scales much larger than the acoustic horizon (lot <C 1 in eq. 
|26| ), the sum of the photospheric terms in equation ([^) (the terms evaluated 
at recombination) is |(/>. Thus, on angular scales much larger than 1° (roughly 
the size of the acoustic horizon) , the CMB anisotropy is a direct measure of the 
gravitational potential on the photosphere at recombination. 

4.1. Angular Power Spectrum 

The angular power spectrum gives the mean squared amplitude of the CMB 
anisotropy per spherical harmonic component. Thus, we expand the anisotropy 
in spherical harmonic functions of the photon direction: 

A(n) = J2 a imYim(n) . (40) 

l,m 

Observations of our universe give definite numbers for the a/ m (with experimental 
error bars, of course). Theoretically, however, we can only predict the probability 
distribution of the a\ m . For statistically isotropic fluctuations (i.e. having no 
preferred direction a priori), the a\ m are random variables with covariance 

{aimdi'm 1 ) = Ci8iii8 mm i (41) 

where the Kronecker delta's make the a; m 's uncorrected (du> = 1 if I = I' and 
otherwise). The variance of each harmonic is given by the angular power 
spectrum Cf, rotational symmetry ensures that, theoretically, it is independent 
of m. 

To calculate the angular power spectrum, we expand 4>(x., r) (and the other 
variables) in plane waves (assuming that the background is flat, K = 0), 

/~ oo 
d?ke- i ^(i>(k,T)= d 3 k^2(-i) l (^ + m(kx)Pi(^(k,r) (42) 
J 1=0 
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where fi = —k ■ n/k with a minus sign because — n is the radial unit vector at 
the origin. Note that many cosmologists insert a factor (2-7r) -3 in the Fourier 
integral going from &-space to x-space; take heed of the conventions when they 
matter. 

In equation (|42] ) we have used the spherical wave expansion of a plane wave 
in terms of spherical Bessel functions ji(x) and Legendre polynomials Pi{x). 
This gives 

„ CO 

A(n) = / d 3 k XX-*)'(2J + l)i^)A,(fc) (43) 
J 1=0 

where, in the instantaneous recombination approximation, 



— 5j(k, T rcc ) + 



jl(kXvec) + kW r (k, T rec )jl(kXr 



ret' 



+ / d X 2d T( f>(k,To-x)ji(kx) ■ (44) 
Jo 

The temperature expansion coefficient takes a simple form in terms of A^: 

aim = H)'4tt / d'kY^Mk) (45) 



where A: is a unit vector in the direction of k. 

To get the angular power spectrum now we must relate A/ (A;) to the ini- 
tial random field of potential or entropy fluctuations that induced the CMB 
anisotropy. Let us assume we have adiabatic fluctuations, for which <p(k,T — > 
0) = A(k). We define the CMB transfer function 

Dm , *S® , (46) 
A(k) 

which depends only on the magnitude of k because the equations of motion have 
no preferred direction. To compute ensemble averages of products, we will need 
the two-point function for a statistically homogeneous and isotropic random 
field, whose variance defines the power spectrum: 

(A(k )A*(k')) = P A (k)5l(k - k>) . (47) 

Note that many cosmologists use an alternative definition of the power spectrum 
with a factor (2-7r) 3 inserted on the right-hand side. The power spectrum so 
defined is greater by a factor (2-7r) 3 than the one given here. [The extra factors of 
(2-7r) 3 in eq. ^2] are then important.] The definition presented here is conventional 
in field theory and it ensures that PA(k)d 3 k is the contribution per unit volume 
in A>space to the variance of A(x). Because of this interpretation, the power 
spectrum is also known as the spectral density. See |]l5| for further discussion. 

For scale-invariant fluctuations (equal power on all scales, as predicted by 
the simplest inflationary models), Pa oc k~ 3 . 
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Combining equations (fll]), (15), and fl4"6|) gives the formal expression for 
the angular power spectrum as an integral over the three-dimensional power 
spectrum: 

Q = 4vr J d 3 kP A (k)D?{k) . (48) 

It is difficult to get a simple approximation to D[(k) that makes this analytically 
tractable, except in the limit of large angular scales where the intrinsic and grav- 
itational contributions to A(n) dominate, with A ~ |0. Then D\ = ^ji(kxrec) 
and the integral may be evaluated analytically for -Pa(^) °c k n ~ A with fixed n. 
When n = 1 (the scale-invariant spectrum), the result gives l(l + l)Ci = constant 
or equal power on all angular scales. The phenomenon of acoustic peaks in 
/(/ + 1)C; is due to the acoustic oscillations in (f> and W 7 which modify the 
potentials at recombination from their scale-invariant primeval forms. 

This presentation has been based on the traditional Fourier space represen- 
tation of the potential and other fields. The physical interpretation of the CMB 
two-point function is somewhat clearer in real (position) space, although a more 
detailed analysis based on Green's functions is needed before one can reap the 



rewards. For an introduction to the position-space perspective see [16|. 



5. Structure Formation in ACDM Models 

This section summarizes briefly the steps needed to go from fluctuation evolution 
in the early universe to predictions of galaxy formation and large scale structure. 
The presentation will focus on techniques rather than model results. However, 
in order to simplify the discussion, we will focus on the favored cold dark matter 
family of models (with curvature and vacuum energy as options) with adiabatic 
density fluctuations from inflation. 

5.1. CDM Linear Transfer Function 

During the radiation-dominated era, CDM density fluctuations grow only loga- 



rithmically. Once the universe becomes matter-dominated at r ~ r e (eq. 24), 
equation (|l6| ) shows that the potential becomes frozen in as long as vacuum 
energy (or other sources of pressure) and curvature are dynamically unimpor- 
tant. Thus, from equation (|14|), the CDM density fluctuation amplitude grows 
linearly with a since a eq : 

6 m (k,r > r e ) = ~jf}_ = -|^(A;r e ) 2 0(fc,r e ) . (49) 

We have used subscript "m" to include all cold matter including baryons, which 
cluster with the CDM after recombination. We have also included the velocity 
potential term 3(a/a)W in equation ( |14| ) as part of S m as it is in Bardeen's 
gauge-invariant density perturbation variable e ||. (The velocity potential term 
is anyway negligible on scales much smaller than the Hubble distance H .) 

To fully specify the CDM density fluctuation spectrum we now need to 
know the gravitational potential, or equivalently the density fluctuations, at the 
end of the radiation era. There is a crucial physical scale, the Jeans length for 
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the radiation gas at r = r e 



UMpc 



(50) 



(Including the effects of baryons decreases this by up to 15% from what is given 
here.) Perturbations of wavelength longer than 2ir/k e have always been Jeans 
unstable, so that <j)(k,T e ) ~ A(k) for k <C k e . Note that only the acoustic 
horizon — not the causal horizon or Hubble scale — enters as a length scale in the 
CDM transfer function. 

To get the amplitude for shorter-wavelength perturbations which have suf- 
fered acoustic damping (i.e. they were Jeans stable for a period of time), it is 
not enough to use the radiation-era solution for the gravitational potential in 
equation (|4^), because the CDM density fluctuations dominate over those of 
radiation before the universe becomes matter-dominated. Instead, we turn to 
equations (|14|) and ( |26| ) for ujt 3> 1 in the radiation era. For the adiabatic mode 
this gives 

3 9^ i— 

S c = ~6 r — 7] rs -A(k ) cos(a;r) — r\ , where uj = kj v3 



-9A(k] 



ln(wr) + C 



Ci(wr) 



r < r e and ujt 3> 1. (51) 



The gravitational impulse created by the acoustic oscillations of the radiation 
fluid causes the CDM to oscillate slightly but as ujt grows the cosine integral 
quickly dies and the CDM grows logarithmically. This growth is important for 
structure formation as it allows the CDM perturbations to dominate over the 
radiation perturbations enough for galaxies to form in less than a factor 1000 
expansion since recombination when the large-scale radiation fluctuations were 
only a few parts in 10 5 . 

Combining results, we obtain an approximate solution for the CDM transfer 
function T m (k) giving the amplitude of CDM density fluctuations relative to the 
primeval potential A{k ): 



T m {k) 



-2(k/k e ) 2 , 
-9[ln(k/k e )+C] 



k <t:k e 

/c » k e 



where 



,r > r e 



o(r) 



<-eq 



A{k)T m {k) 



(52) 



(53) 



Bardeen et al. [17| give a fitting formula for T m (k)/k 2 , which, when normalized 
to 1 as k — > 0, is conventionally called the CDM transfer function. However, 
given the measurements of the gravitational potential made possible from the 
CMB, it is preferable to retain the normalization of equation fl52| ) with the 

potential A(k) factored out in equation (|53|). 

The CDM linear power spectrum follows by analogy with equation ([47]) for 
the potential. The result, for a matter-dominated universe, is 



Pin{k, t) 



o(r) 



Tl{k)P A {k) . 



(54) 



18 



Edmund Bertschinger 



CO 

.£= 

O 
Q_ 



CD 
O 




-1.5 -1 
log 10 k (h/Mpc) 



Figure 1. Theoretical CDM power spectra compared with the decon- 
volved power spectrum of galaxies from the Las Campanas Redshift 
Survey [19]. The error bars are strongly correlated at large scales. 



The constant-curvature scale-invariant spectrum of potential fluctuations Pa oc 
k~ 3 was first proposed by Harrison IS] for its naturalness properties long before 



inflationary cosmology gave it a concrete basis. With this spectrum, the CDM 
density spectrum has the well-known limiting behavior 

P ™W " {k 2 *Mk/ke)? , ^ 

The power spectrum is useful for calculating the variance of the smoothed 
density field. If we convolve 5(x, r) with a spherical window function whose 
Fourier transform is W(k) normalized so that W(0) = 1, then 



J d 3 kP(k)W 2 (k) . (56) 



From this, we see that the CDM density fluctuations have logarithmically diver- 
gent power at small scales. Thus, structure formation in a CDM model proceeds 
by hierarchical or bottom-up clustering with small scales collapsing by gravita- 
tional instability before being incorporated into larger objects. However, the 
slowness of the logarithm means that the structure formation proceeds rapidly 
from small to large scales. 

Figure |l] compares the linear power spectra extrapolated to r = To (i.e. 
redshift 0) for several flat CDM models. The curves are all normalized so that 
the density smoothed with a 8 h~ l Mpc spherical "tophat" has standard devi- 
ation erg = 0.9 in agreement with observed galaxy clustering for a "linear bias 
factor" b = cig 1 = 1.1. The ACDM model with (Q b = 0.05, (l c = 0.30, ft A = 
0.65, h = 0.65) has the most power on large scales, possibly too much (although 
an acceptable fit can be found with slightly higher f2 c ). The SCDM model with 
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(f2b = 0.05, O c = 0.95, h = 0.50) has too little power on large scales. When 
combined with the CMB anisotropy, this provides a strong argument for dark 
energy (e.g. a cosmological constant) |p(3| . 

The curves labeled "BBKS" show the linear power spectra for the ACDM 
and SCDM models in the limit 0,f, = (with f2 c = 0.35 and 1.0 for the two mod- 
els, respectively). They are obtained from the fitting formula given by equation 



(G3) of |Tj and agree well with full calculations using the COSMICS code f21fl 
which was used for the other curves in Figure [l|. Baryons reduce the amount of 
small-scale power relative to scales k < k e because they cannot experience grav- 
itational instability on scales smaller than the photon-baryon Jeans length until 
after recombination. The oscillations in the ACDM power spectrum are also 
due to baryons and are similar to the cosine integral contributions to 5 c {k) in 



equation (51). They arise from the acoustic oscillations of the baryons between 
r e and r rec . When the power spectra are normalized at small scales (as they 
effectively are with ag), increasing the baryon fraction increases the amount of 
large-scale power. 

Aside from the baryon effect, the main parameter affecting the CDM spec- 
trum (and the only one in the BBKS formula) is k e . As may be checked using 
equation (|50|), the peak of the CDM power spectrum is closely equal to k e . 
Measurements of galaxy clustering actually constrain k e /h = 0.1T Mpc -1 where 
r = Vt m h is a commonly used parameter. However, as Figure |l] shows, a sec- 
ond parameter, the baryon fraction £lb/£l m , also significantly affects the CDM 
spectrum. 

5.2. Evolution of Large Scale Structure 

The results given above concerned the shape of the CDM power spectrum, which 
is frozen in after the universe becomes matter-dominated. However, the ampli- 
tude grows as a(r) as assumed in equation (|53|) only if Q m = 1. While this is a 
good approximation soon after recombination, it breaks down at small redshift 
when Q m < 1 today. In this case we must include dark energy, curvature, and 
any other important contributor to the expansion rate of the universe. 

The evolution of CDM density fluctuations for r 3> r e is given by equations 
( |l4| ) and ( |2(i| ) ignoring the radiation terms. On subhorizon scales for (j) 2 <C 5^ <C 
1 we have 

(X 3 $ 

dr8 m + -d T 5 m = i7rGa 2 pmS m = Ti^mHl — . (57) 
a 2 a 

Because the coefficients are spatially homogeneous, the time and space depen- 
dences of the solution factor: <5 m (x, r) = 5 + (x)D + (t) + <5_ (x )L>_ (r). For 
O m > 0, D + grows with time and L>_ decays. For Q m = 1 the growing solution 
is D + = a. For £l m < 1 the solutions depend on the curvature K (sometimes 
expressed as VLk = —K/Hq) and on the equation of state of the dark energy. 



See Peebles |2j and Peacock |ll| for examples. 

Peculiar velocities offer another probe of large-scale structure that is sensi- 
tive directly to all matter, dark or luminous. The divergence of peculiar velocity 
v = dx/dr depends on the growth rate of density fluctuations: 

9 m = V • v m = -d T 5 m = - C^^f) aH5 m for 5 2 m « 1 . (58) 
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Measurements of galaxy velocities actually give v/aH, so that the ratio of ve- 
locity divergence to density fluctuation gives the logarithmic growth rate of 
density fluctuations. An empirical fit to a range of models gives din D + / din a = 
f{Q m ) ~ f^ ' 6 - In practice, the galaxy number density may be a biased tracer 
of the total matter density, conventionally parameterized by the bias factor b: 
^galaxies = bS m . Comparing galaxy clustering and velocities then allows determi- 
nation of WJ/b. 

5.3. Numerical Simulation Methods 

Linear theory breaks down once dark matter clumps become significantly denser 
than the cosmic mean. Once this happens, full numerical simulations are needed 
to accurately follow the nonlinear dynamics of structure formation. Moreover, 
radiative processes, star formation, and supernova feedback have a strong effect 
on the baryonic component once galaxies begin forming. Although these pro- 
cesses can only be included in simulations in a phenomenological way, numerical 
simulations are still the best way to advance our theoretical understanding of 
structure formation deep into the nonlinear regime. 

Here we present only a brief synopsis of numerical simulation methods. For 
a comprehensive review of both the methods and their applications see |^3| . 

Particle methods are used for representing the dark matter. The phase 
space of dark matter is represented by a sample of particles evolving under their 
mutual gravitational field: 

dxi _ dvi 2 2 

Vi , —r- = Vi — Vcp , V = Air&a 



dr dr a 



^mjd-D{x - Xi) - p m {r) 

(59) 

where <5d is the Dirac delta function. In practice, the particles are spread over 
a softening distance to avoid unphysical two-body relaxation due to large-angle 
scattering. (A simulation may use particles of mass 10 9 M , leading to un- 
physical tight binaries if the forces are not softened.) The art of dark matter 
simulation is mainly in calculating forces much faster than 0(N 2 ) for a soften- 
ing distance several orders of magnitude smaller than the simulation volume. A 
variety of algorithms and codes are available; see [23] for details. 



The initial conditions for dark matter positions and velocities are given from 
linear theory by the Zel'dovich approximation J7||, which relates the comoving 
position to its initial value q, a Lagrangian variable (i.e. one that labels each 
particle) : 

x(q, r) = q + D + (r)$(q) , V 3 • f = -S m . (60) 

The initial density fluctuation field 5 m (x, r) is obtained as a sample of a gaussian 
random field with appropriate power spectrum. 

Baryon physics is added by numerical integration of the fluid equations (|l7|) 
and (|i"8|) with the addition of heating and cooling, ionization, etc. There are 
two classes of methods in widespread use in cosmology: smoothed-particle hy- 
drodynamics (SPH) and grid-based methods. SPH is an extension of the N-body 
dynamics of equations ( |59| ) to include pressure forces and thermodynamics. It 
has the advantage of concentrating the spatial resolution and computational 
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effort where the baryons are densest, with the drawback of relatively large nu- 
merical viscosity and diffusivity owing to the particle discreteness. Grid-based 
methods use finite-difference techniques to approximate the fluid equations on 
a (usually) regular lattice. High accuracy methods developed for aerodynamic 
and other applications provide excellent resolution of shock waves, but with a 
resolution limited by the grid. Adaptive mesh refinement techniques are now 
being used to refine the resolution where needed, resulting in the ability to follow 
structure formation to much higher densities than previously possible j24|. 

We conclude this introduction to numerical simulation methods with a mo- 
saic of gravitational potential and dark matter density from a high-resolution 
simulation of the mixed cold+hot dark matter model. Figure [2] shows the linear 
theory and nonlinear simulation results. On the scales shown, ACDM would be 
almost indistinguishable. 

The gravitational potential shown in Figure § is a statistical fractal at the 
end of inflation, viz. a gaussian random field with power spectrum Pa oc fc~ 3 . 
During the radiation era, acoustic damping smooths the potential on scales 
smaller than k~ l . Because the density fluctuations are related to the potential 
by two spatial derivatives, the corresponding linear density field is dominated by 
high spatial frequency components. Indeed, if not for the smoothing imposed by 
finite spatial resolution, the variance of density fluctuations would be logarith- 
mically divergent. It is instructive to compare the post-recombination density 
with the inflationary potential. The density is almost proportional to the po- 
tential but with opposite sign. This is because the density/potential transfer 
function (|52|) varies only logarithmically on small scales. 

The nonlinear evolution of the density field causes the mass to cluster 
strongly into wispy filaments and dense clumps; most of the volume has low 
density. The intermittency of structure is a sign of strong deviations from gaus- 
sian statistics. 



5.4. Press- Schechter Theory 

Numerical simulations show that most of the mass falls into dense self-gravitating 
clumps. As a crude approximation, let us suppose that the mass profile of every 
clump is spherically symmetric. In this case every gravitationally bound shell 
will eventually fall into the center of the clump, at a time that is dependent 
on the enclosed mass and energy of the shell. For linear growing mode initial 
conditions, collapse occurs when linear theory would predict that the mean 
density contrast interior to the shell is 5 C = 1.686 [[D|. Thus, in the spherical 
model, the distribution of clump masses is simply dependent on the statistics of 
the smoothed linear density field. 



Press and Schechter [26] devised an elegant method for computing the mass 
distribution of gravitationally bound clumps using the simple spherical model. 
Their procedure is first to smooth the linear density field by convolution with a 
spherical "tophat" of radius R enclosing mean mass M(R) = (47r/3)/9 m i? 3 : 

d 3 x'W R (\x-x'\)8 m (x), W*(0 = J^{J; r r <*> (61) 

According to the Press-Schechter ansatz, the regions of the initial density field 
with 5 > 5 C = 1.686 at time r have collapsed into objects of mass at least M. It 
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Figure 2. Potential and density for a mixed cold+hot dark matter 
model with £l v = 0.2 ||25|]. A slice of size 50/i _1 Mpc is shown. Upper 
left: potential at the end of inflation. Upper right: potential at the 
end of recombination. Lower left: density at the end of recombination. 
Lower right: density at redshift 0. 
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is easy to see that this would be true for a single spherical density perturbation; 

the genius is applying this to the entire density field. 

Under this model, the fraction of mass in the universe contained in clumps 

more massive than M is 

1 f°° dn 
— \ dMM— = P(S > S C \M, t) . (62) 
p m Jm dM 

On the left-hand side, {dn/ dM)dM is the number density of clumps of mass in 
(M, M + dM) ; the integral gives the total mass density in clumps more massive 
than M. Dividing by p gives the mass fraction. On the right-hand side, P(5 > 
6 C \M, t) is the probability that a randomly chosen point in the linear density 
field has 5 > 6 C at time r, where M defines the smoothing scale of 5. For a 
gaussian random field of density fluctuations, this probability is simply 

P = Lrfc (63) 



2 \<t>J2, 

where erfc is the complementary error function and a 2 = a 2 (M,t) = (5 2 ) is the 
variance of the smoothed density. From equation (|62|) we can deduce dn/dM by 
differentiating with respect to M. 

Now, as M — > 0, a{M) grows without bound for CDM models with loga- 
rithmically diverging power on small scales, so that P{S > S C \M = 0,r) = i. 
Half of the mass initially is in over dense regions, half is in under dense regions 
5 < which in linear theory can never collapse. However, this is plainly unphys- 
ical, as simulations show that virtually all of the mass, including that in initially 
underdense regions, accretes onto dense clumps. Thus, Press and Schechter sug- 
gested doubling dn/dM so that all of the mass resides in dense clumps. Thus, 
their formula for the mass distribution of clumps is 

dn _2p dP 

dM ~ MdM ' ^ ' 

Although realistic gravitational collapse is far more complex than imagined 
in this derivation, equation (m) gives surprisingly good agreement with clump 
masses measured in cosmological N-body simulations j27j. The Press-Schechter 
model and its generalizations have proved to be a very useful tool for comparing 
ab initio theories (power spectrum and cosmological model) with observations 
of galaxies and galaxy clusters. 



6. Conclusions 

These lecture notes have presented a few of the theoretical elements useful for un- 
derstanding the microwave background and structure formation. Cosmology is 
rapidly changing with the advent of precision measurements of CMB anisotropy, 
with the observation of galaxies at z > 3, and with large new redshift surveys. 
Theoretical research has become more phenomenological, with a focus on pro- 
viding the framework for interpreting and analyzing the new data. Despite 
this trend — or perhaps because of the high quality of new data — theoretical 
cosmology still requires further development of the basic processes of structure 
formation. Either way, I hope that the material presented in these lectures will 
be useful to cosmology graduate students. 
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